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Abstract 

In this work we study a microtubule (MT) model, whose length is regulated by 
the action of processive kinesin motors. We treat the case of infinite processivity, 
i.e. particle exchange in the bulk is neglected. The exact results can be obtained for 
model parameters which correspond to a finite length of the MT. In contrast to the 
model with particle exchange we find that the lengths of the MT are exponentially 
distributed in this parameter regime. The remaining parameter space of the model, 
which corresponds to diverging MT lengths, is analyzed by means of extensive Monte 
Carlo simulations and a macroscopic approach. For divergent MTs we find a complex 
structure of the phase diagram in terms of shapes of the density profile. 


1 Introduction 

The cytoskeleton is a dynamic network of biopolymers which determines the shape of en- 
karyotic cells and serves at the same time as a transport network. One distingnishes three 
different types of cytoskeletal filaments (actin, microtnbnles and intermediate filaments). Mi- 
crotnbnles (MTs) consist of a and (3 tnbnlin snbnnits and form cylindrical strnctnres with 
high bending rigidity |T]. The dynamic properties of MTs allow them to take on arbitrary 
length-distribntions, which are adapted to the different cell types and shapes. Therefore 
it is of fnndamental interest to identify the mechanisms, which control the length distribn- 
tions of MTs. These regnlatory mechanisms are carried ont by polymerases, which snpport 
MT growth, or by depolymerases, which are able to induce fast depolymerization events, 
so-called catastrophes [2j. The stabilization of growing MTs has strong influence on their 
length distribution, in particular under spatial conhnement [3]. Dynamically stabilized MTs 
are able to connect the cell center and membrane without losing their ability to adapt to 
different cell shapes. Other proteins have been identihed which are able to enhance the 
depolymerization of MTs. This is the case for Kip3 a member of the kinesin-8 family with 
extremely high processivity. Kip3 moves to the plus end of MTs where it is able to enhance 
the depolymerization activity. This process is of great importance in order to position the 
mitotic spindle in yeast [H El E]. 
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The length regulation mechanism induced by Kip3 has been addressed theoretically in 
[3 El El [10], by using one-dimensional driven diffusive processes. These models capture the 
stepwise directed motion of the kinesins and describe their interaction by mutual exclusion. 

From a theoretical point of view driven diffusive processes are stochastic models of in¬ 
teracting particles, which have been studied typically on lattices with a hxed system length 
(i.e. with a hxed number of sites) or on Z [TT]. The totally asymmetric simple exclusion 
processes (TASEP) is the archetype of this kind of models. One of well-studied cases is the 
TASEP with open boundaries, i.e. on a hnite chain, particles are injected at the one end 
and extracted at the other end of the chain [I211I3]. We refer to this realization of the model 
as open TASEP. Recently the TASEP or its related processes with varying system size have 
been intensively studied, motivated by other biological applications or modeling queues with 
exclusion principle [El [El lEl lEl lEl lEl [ 2 Q| as well as the MT models [71|8l|9l [El El] . 

In this work we study a variant of the TASEP, which was introduced in [5] as a “sim- 
plihed” model of length regulation via kinesin motors. The parameters characterizing the 
model are the motor protein input rate, the polymerization rate and the depolymerization 
rate. Exchange of motors in the bulk, i.e. Langmuir kinetics, is not considered, which 
seems to be justihed in view of the high processivity of Kip3. (Langmuir kinetics was also 
introduced in the full model laio 

We denote the length of the MT by L, where each site is labeled by an integer j from the 
minus end j = 1 to the plus end j = L, see Fig. (a). Molecular motors occupy a single site 
and interact via mutual exclusion. We write Tj = 1, 0, if the site j is occupied or unoccupied, 
respectively. Therefore the state space S consists of configurations r = > 0), 

where we write r = 0 for L = 0. A motor at site j walks to site j -|- 1 with rate 1, if site 
j -|- 1 is not occupied by another motor. A new motor enters the system with rate A (in [H], 
this is denoted by p_) at the minus end if the hrst site exists and is not occupied. 

In contrast to the open TASEP the lattice size which represents the length of the MT, 
is not hxed. At the plus end of the MT, polymerization occurs, i.e. a new site is attached 
to the lattice with rate 7 . If the MT is completely depolymerized, i.e. L = 0, a new site is 
created with rate 7 . Depolymerization occurs with rate 6, if the site L is occupied. In this 
case the last site is removed together with the particle that occupies this site. Note that we 
do not take into account the cell size, i.e. no upper bound of the MT length L is introduced. 

In this work we present an exact analytical approach to the stationary properties of the 
model for hnite MT lengths. We also use a phenomenological approach and Monte Carlo 
simulations for the case of diverging MTs. We show in section (and in for details) 
that the stationary measure can be written in a matrix product form, which is a “grand 
canonical” version of the open TASEP. By evaluating the normalizability of the stationary 
measure, we determine the area in the parameter space where the MT length converges. In 
this case, the average MT length can be exactly calculated by using some known results 
from the open TASEP. In order to analyze the divergent case in section we give some 
analytical arguments: For divergent MTs the macroscopic density prohles can be of shock 
or rarefaction type, depending on the parameters of the model. We also check by Monte 
Carlo simulations that the predicted density prohles are actually realized. In Section we 
summarize this work and give some concluding remarks. 
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Figure 1: (a) Schematic illustration of the model dynamics: Each binding site of the MT, 
which is represented as a lattice site, is occupied by a single motor or empty. New motors 
are injected only at the — end, which can be regarded as reservoir with density A. At the 
+ end, polymerization and depolymerization occurs. We assume that the depolymerization 
is allowed only if the last site is occupied by a motor. (b,c) The phase diagram in terms 
of converging MTs and diverging MTs (c = The convergent phase appears when 

7 < I (b). For 7 > I (c), MTs always diverge. 

2 Stationary distribution 

Depending on the growth rate 7 , one observes a stationary length distribution of the MTs 
(7 < 7 c) or a divergent MT length (7 > 7c). In |9] the critical rate 7 c has been established 
for different input and depolymerization rates 

r A(l-A) (A<|,A<5), 

7 c= <5(l-<5) (5<iA>5), ( 1 ) 

U (A,5>^). 

Figure (^, 0 ) shows the phase diagram in \-6 plane for given 7 , which is equivalent to this 
result. The convergent phase is the region where both S and A are larger than the critical 
value c = (a solution to c(l — c) = 7 ) for 0 < 7 < 

The form of 7c is mathematically the same as the stationary current of the open TASEP 
with fixed lattice size, where particles are injected with rate A and extracted with rate 6 

[T 2 I |13] . In [13] it has been shown that the stationary weights of the open TASEP on static 

lattices can be represented by (in general) inhnite dimensional matrices Xi,Xo and vectors 


(IFI, \ V) that satisfy the following algebra: 

XiXo = Xi+Xo, (2) 

A(IF|Xo = (IF|, (3) 

5X,\V) = \V). (4) 

For a given system length L, the stationary probability of finding a configuration ti - ■ - tl is 
represented as [13] 

P{n---n) = Xw\Xr,---XrJV) (5) 
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with the normalization 


Zi= Y.{W\X^,---X^,\V) = (W\(X, + X„)‘-\V). (6) 

Tj=lfi 

From ([^ it is evident that the stationary current is calculated as 

J= {r,(l-r,+i)) = ^(PF|(Xi + Xor'XiXo(Xi + Xo)^-*-'|l/) = (7) 

In [13] it was also shown that the stationary current approaches the values described in Q 
as L —)■ cx). 

Let us go back to the MT model. We note that, in contrast to the open TASEP, the 
MT length L is also a variable. A “grand canonical ensemble” of the open TASEP gives a 
stationary measure of this process, where the polymerization rate 7 plays the role of fugacity 
[ 22 Q(seegfor details): 

F(0) = 1, F(ri---ri)=7^(hF|X.,---X.jE). (8) 


A similar realization of the grand canonical ensemble of the open TASEP has been pre¬ 
viously found in [23] . Polymerization and depolymerization at the right boundary lead to 
mathematically the same condition as in j23j^ This recipe has been also applied in order to 
construct an exact stationary solution to a queueing process with exclusion [18], where the 
arrival rate plays the role of fugacity. However, constructing a matrix product form with 
some parameter as fugacity is not always possible. One counterexample is given by another 
queuing process with high and low priorities of customers PDl, where it seems difficult to 
hnd a matrix product solution although this model can be regarded as an exclusion process 
with varying length. 

The form ([^ is a direct consequence of the normalizability of the stationary measure. 
In order to check for convergence of the normalization we hrst note that the stationary 
measure for a given length L is given by 


TdS T,=1,0 

(length=L) 


The normalization of the open TASEP was calculated as [T3] 


^L = 


^L!(L-j + l)!l/A-l/5 




[A(l-A)]-‘ (A<l.A<i). 
V (A.i>|). 


( 10 ) 


^Here 7 '^ corresponds to the states where the number of lattice sites is L (not the number of parti¬ 
cles). When the polymerization rate 7 ^ depends on the MT length L, a simple generalization of Q with 
replacement 7 '^ 7071 • • • 7l-i describes the stationary measure. 

^In the model of [23], the left end is also dynamical, where a site occupied by a particle is attached. When 
the leftmost site is empty, it is removed from the system. 
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From this behavior we hnd that 

2 = ^F(r) = 5^7^Zz, (11) 

tGS L>0 

is hnite when 7 < 7 c with equation Q. In other words, for 7 < 7 c the stationary measure 
^ is normalizable and therefore ^F{t) gives the exact stationary distribution. The critical 
value 7 c corresponds the case where limi^oo = 1. 

The grand canonical partition function of the open TASEP is known to have the form 
[251126] 


4X6 


z= _ !l9i 

(^r^4^-l + 2A)(x/r^4^-l + 25) ^ ^ 

for 7 < 7 c. Furthermore, the stationary distribution of the length L is given by Z, 

and its mean value converges to 


(i) = Ei 

L>1 


7''Zl 3 , , 

=7 log2= 


47(x/1 - 47 - 1 + A + 5) 


47(yr^4^ - 1 + 2A)(^r^4^ -1 + 26)' 


(13) 


Except for the cases where A is close to the phase boundary A = c, the average length does 
not take a large value (typically e.g. (L) < 10), as shown in Fig. [^(a). Moreover, the lengths 
L are (except for very small values) exponentially distributed, as one can gather from (10): 


F(L) = ^^Zl/Z 


7 


A(l-A) 


7 


5(1-5) 

47 )^ 


(A: c<A<|,A<5), 
(B: c<(5<|,A>(5), 


(14) 


(C: A>|,5>|). 

One can distinguish three subphases in the convergent phase as shown in Fig. |^(b), according 
to the asymptotic behavior of the length distribution. 


The formula (13) differs from Equation (8) of |9|, where the Langmuir dynamics in the 
bulk is taken into account. These results underscore the importance of Langmuir dynamics 
for length regulation of MTs in inhnite volume. 

The inflow and outflow of lattice sites should be balanced in the stationary distribution. 
In other words, the “tip density” p+ = (r^) must satisfy 7 = 5p+. The matrix product 
solution is compatible with this relation. In fact we have 

P+ = 4 + X„)‘-'X,\V} = I (15) 




L>1 


thanks to relation (|^. 

In this section, the matrix product form enabled us to show the condition 7 < 7 c with 
([^ for the convergence of the MT length, and to calculate the length distribution in the 
convergent phase. However, it does not give a rigorous result for 7 > 7 c, since in this case 
no stationary state is reached and the length of the MT diverges. This implies for example 
that (15) is no longer the correct tip density when the MT length diverges. Thus we shall 


rely on Monte Carlo simulations and a phenomenological approach to hnd density prohles 
in the divergent phase in the next section. 
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Figure 2: (a) Average MT length (L) ( [13) ) vs A in the stationary distribution. The values 
of S are indicated in the panel, and 7 = 0.16. The average length diverges as A approaches 
the critical value c = _ q 2 . (b) Subphases of the convergent phase. In the three 

regions A, B and C, the length distribution exhibits the different asymptotic behaviors as 
Equation (14). 


3 Divergent phase 

In this section we will characterize the states in the divergent phase on a macroscopic level, 
by using the conservation of mass. 

Macroscopically we can identify domains of different density p{x), in which the current 
J{x) is given by J{x) = p{x){l — p{x)). In general the current in different domains may 
differ, which means that mass is transferred from one domain to another. Let us consider 
the typical two waves appearing in the TASEP [21]. If mass is transferred from a high to a 
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Figure 3: Macroscopic density prohles for particles that move in positive x-direction. (a) 
Particles move from a high density domain (pr) toward a low density domain {pe)- A rarefac¬ 
tion wave evolves. The left (right) end of the slope part is moving with velocity ve = 1 — 2pi 
{vr = 1 — 2pr). (b) A shock wave is established if particles move from a low toward a high 
density domain. 
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low density domain, this leads to rarefaction density profiles 


[pi {x <1- 2pi) 

P(^) = S {I-2pi < X <l-2pr) 

Pr (x > 1 2p^) 

with Pi > pr, see Fig. (a). Shock wave density prohles 

Pix) = I '’■ = >’'■> 

\ p, (x > II,) 


( 16 ) 


(17) 


are observed if pf < p,, see Fig.^(b). The rarefaction and shock wave positions are stationary 
in the frame in which the spatial position j is scaled by the time t as x = j/t. 

In earlier works [19] and [20], density prohles of the qneneing processes with exclnsion 
were analyzed. It tnrned ont that the macroscopic density prohles can be rarefaction or 
shock type, respectively. We shall show that both types are also realized in the MT model, 
depending on the parameter valnes. 

We are interested in the macroscopic density prohle of the MT p{x) (0 < x < n+), where 
denotes the velocity of the MT’s pins end. For onr analysis it is important to notice 
that the actnal density p_|_ at the pins end is, in general, diherent from the density R of the 
macroscopic density prohle at the pins end of the MT: 


R := lim p{x). (18) 

X 

The tip velocity and density are related to each other via 


v+ = 'y- 6p+ (19) 

(c/. Eqnation (1) of [9]), since the ehective depolymerization rate of the MT is given hj 6p+, 
which is identical to particle cnrrent in the frame of the tip, i.e. in the co-moving frame. 
Therefore, the particle conservation law in the co-moving frame is expressed a^ 

R(l-R)-R(j-Sp^) = 5p+, (20) 


where the hrst term denotes the cnrrent generated by particle hopping, the second term is 
the cnrrent indnced by moving the frame. In [9], snbphases in the divergent phase, i.e. EX, 
MC and IN phases, were fonnd according to the form of p+. Here we investigate which kind 
of macroscopic density prohles can be observed in the snbphases, and we inversely determine 
the phase bonndaries by considering conditions for existence of the density prohles. 


3.1 EX phase 

Let ns begin with the assnmption that the tip density p+ prefers to become identical to the 
right end of the macroscopic density R. In this case we have 

Pv = l-^ (21) 

^This is equivalent to Equation (2) in |9], where R is denoted by pb- 
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due to the particle conservation (20) with = R. This case corresponds to the EX phase 
as labelled in j9]. The macroscopic density prohle is not always simply flat, which was not 
emphasized previously. It can become a rarefaction wave (16) or a shock wave ( [ItI ) with left 
and right domain densities as = A and pr = p+- In the EX phase there are in general four 
types of macroscopic density profiles: 


A (0<x<l-2A), 


d(^)|EX-I= ^ 

(l-2A<x< l-2p+), 

[ P+ (1 - 2p+ < X < 1;+), 

(22) 

/'(^)Iex-ii = ] 

r (o<x<i-2p+), 

[ P+ (1 - 2p+ < X < n+). 

(23) 

EX-III = P+ (0 < X < n+). 

(24) 

/'(^)Iex-iv = ] 

f A (0 < X < n<j) 

[ p+ (n, < X < n+). 

(25) 

with the shock and tip velocities 
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A, v+-'j 5(1 

(26) 


Because of the definition of the EX phase (21), the boundaries between this subphase and the 
other two subphases can also be given analytically by simple considerations: The transitions 
from the EX phase to the others occur when the relative size of the density domain p{x) = p+ 
shrinks i.e. 1 — 2p_|_ < n_|_ (in the case of the rarefaction wave) or Vg < n+ (in the case of the 
shock wave). We obtain 5 = 1 — ^ and X = S, respectively. 

In order to determine the substructure of the EX phase, one has to check where each 
form (I,..., IV) of the density profile is realized. This can be done by evaluating p^ ^ Pr, 
Us ^ 0, 1 — 2p+ ^ v+ and 1 — 2A ^ 0. The result is summarized in Fig. Simulation results 
for the four sub-subphases are shown in Fig. (a,b). 


3.2 MC phase 

As we discussed in the previous subsection, a phase transition from the EX phase occurs, 
i.e. the right plateau part pr = p+ (with Equation (21) of the rarefaction wave (A > 1 — 
vanishes, when 6 becomes larger than 1 — y/j, see Fig. Then the plus end breaks into 
the slope part of the rarefaction wave, and therefore the macroscopic density profile has the 
following structure: 


p{x)\ 

p{x)\ 


MC-I 


MC-II 


A 

l—x 
2 

1 — X 


(0 <x < 1-2A), 
(1 — 2A < X < n+), 

(0 < X < n+). 


(27) 

(28) 


The phase boundary between the two sub-subphases is given by 1 — 2A = 0, i.e. A = ^. Let 
US assume the tip density has a stationary value p+. This true tip density is not the limit of 
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Figure 4: Phase diagrams including subphases of the divergent phase. The structure of the 
phase diagram depends on the value of the growth rate 7 , i.e. (a) 0 < 7 < |, (b) | < 7 < |, 
(c) ^ < 7 < 1. The phase boundaries c, d are given by c = and d = 1 — The 

boundaries between EX-I and EX-IV and between EX-III and EX-IV are given as A = 1 — 
and A = respectively. 


the macroscopic density profile (18), but related as 

R = p{v+) = VV- 


(29) 


Here is the tip velocity which is written in terms of as (19). Now we have obtained 
one relation between and R: 

1 - (7 - 6p+) 


R = 


(30) 


In order to calculate p+ we need one more relation, which we can obtain from mass conser- 

(31) 


vation, i.e. Equation (20). Solving them, we find 

(1 - 


p+ = 


6 


v+ = 2y/y-l. 


Simulation results for the two sub-subphases are shown in Fig. |(c,d). 

When (31) becomes less than 1 — 2A (i.e. A < 1 — a/ 7 ), the slope part vanishes and 
we enter the last subphase. 


3.3 IN phase 

As we discussed in the previous subsections, a phase transition from the EX phase occurs 
when 6 becomes larger than A, and another one from the MC phase occurs when A becomes 
less than 1 — ^ 7 , see Fig. In the last subphase, the right plateau pr = p+ (21) of the 
shock prohle (or the slope part of the rarefaction profile) vanishes, and only the left plateau 
remains, which was denoted as IN phase in [9]. The macroscopic density profile is flat in the 
entire region as 


p{x)\m = X {0<x<v+), 


(32) 
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Figure 5: The macroscopic density profiles p{x) in (a) EX, (c) MC and IN phases. We 
compare the results of Monte Carlo simulations (dots) with the predictions by the hydrody¬ 
namic approach (lines). The microscopic density prohles pi in the vicinity of the plus end 
{i = L,L — 1,...) are also shown for the (b) EX, (d) MC and IN phases. We observe no 
boundary layer in the EX phase (b), but the true tip density pi = p+ (indicated by arrows) 
is different from the macroscopic density R (indicated by lines) in the MC and IN phases (d) 
[9], see Equations 


(33) and (34). The chosen parameter values (A, 7 , 5) are indicated in 
the panels (a) and (c), and we used the same values and markers for (b) and (d), respectively. 
As the initial condition we chose L = 0. The results shown in (a) and (c) were obtained 
by averaging 10 000 snapshots at time t = 15 000, and for (b) and (d) we averaged over the 
time window 10 000 < f < 15 000 of 10 000 simulation samples. 
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see Fig. (c). The tip density p+ and the velocity are given as 


_ A(l-A-7) 


= 7- 


A(1-A-7) 

r^A 


( 33 ) 


from (20) with 


R = X, 


(34) 


see Fig. |(d). 

4 Discussions 

In this work we stndied a model for length regnlation of MTs by molecnlar motors introdnced 
in Pi- We restricted onr discnssion of the model to the case where the particle conservation 
is satished in the bnlk, i.e. attachment and detachment of motors in the bnlk are not 
considered. This limit of the model has been snggested in [7] and it has been argned that it 
is relevant for the length regnlation by kinesin 5 dne its high processivity. 

We gave an exact stationary solntion of this model in the parameter regime where the 
average length of the hlaments is hnite. In this regime, we observe exponentially distribnted 
hlament lengths, similar to nnregnlated hlaments. The typical lengths of the MTs in this 
regime are rather short, except for parameter valnes that are at the limit to the divergent 
case. This resnlt implies that length regnlation of the hlament is either generated by the 
bnlk exchange of molecnlar motors, or by varying the inhnx of motors in the vicinity of the 
transition from convergent to the divergent state. 

We also investigated the divergent phase by using a phenomenological approach and 
Monte Carlo simulations. The phenomenological approach enabled us to determine ana¬ 
lytically a rich phase diagram. For inhnitely growing MTs one can distinguish between 
subphases where the entry (IN) and the exit (EX) induce macroscopic domains of constant 
densities. Additionally one observes a phase (MC) where the tip is connected to the slope 
part of a rarefaction wave. The distinction of these three phases had already been made 
in |9]. But in this work we explored the complete substructure of the divergent phase and 
possible shapes of density prohles by hydrodynamic approach, and found good agreement 
with Monte Carlo simulations. 

In view of our results it is of great interest to discuss the biological relevance of different 
regulatory mechanism of the system. In [2] it was argued that processive motors are used for 
MT length regulation since they establish a direct link from the cell center to the periphery. 
For free hlaments, however, our results indicate that length regulation is more robust, if 
the run length of molecular motors is hnite. Therefore it would also be of great interest to 
see under which conditions the positioning of the microtubule organizing center would be 
possible in a hnite volume, as it was observed for microtubules that are pulled by dynein 
motors attached to the boundaries in artihcial cells [28] . 
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A Proof of stationarity 


Here we prove that (|^ gives a stationary measure of the model. The proof is very similar 
to the paper [13], where the matrix product ansatz was firstly introduced in the stochastic 
interacting particle system. To make this work self-contained, we write the proof in detail. 


The master equations concerning r = 0, 0 and 1 are 

|F,(0) =0F,(1) - 7F,(0), (35) 

|F,(0) =0F,(O1) + 7P,(0) - (A + 7)P(0), (36) 

|P(1) =0P(11) + AP(0) - (7 + S)Pt(l). (37) 

For general configurations ti - ■ ■ tl with L >2, the master equation is written as 

|n(Ti • • • Ti) (38) 

= A(2t,- l)P(0T2...Ti) (39) 

L-1 

+ ■ ■ ■ Tj-ilOTj+2 ■■■Tl) (40) 

(1 - TL)'jPt{Ti ■ ■ ■ tl-i) - 'yPtin ■■■Tl) (41) 

-h 5Pt{Ti ■ ■ ■ tlI) - TLdPtin ■ ■ ■ tl-iI). (42) 


It is convenient to clearly distinguish the terms “stationary measure” and “stationary 
distribution” as defined in [27|. A stationary measure is a function F in S' such that the 
right-hand sides of all the master equations become 0 by substitution Pi(r) = F(r). When 
F is normalizable, i.e. Xlrss F(r) is finite, P(r) = F{T)/J2res F{t) gives a stationary 
distribution. On the other hand, when diverges, the stationary measure loses its 

physical meaning m 

Let us find a stationary measure in our case. Without loss of generality we set F(0) = 1. 


We soon notice that F(l) = ? from (35). One wants to know the form of F(0). However we 


cannot find it by solving simultaneously (36) and (37), where apparently we need information 


of F(01) or F(ll). Similarly, to know the form of F{ti .. .tl) of arbitrary configurations 
Tl... Tl, we need information of configurations with length L-|- 1. This fact prevents us from 
systematically finding a stationary measure, and thus an ansatz is required. 

We now prove that the ansatz form actually gives a stationary measure. Checking 


that (36) =0 and (37) = 0 is not difficult. For the general case, the term (39) becomes 


-(l-2r07^(W^l^r,---X,, 


iv 


(43) 


thanks to (1^. Each summand of ( |40[ l becomes 

+ (1 - 2T,)',‘-(W\Xr, ■ ■ ■ Xr^_,Xr,,, ■ ■ ■ X^^V) 
- (1 - 2r^+,)7yW'|A'„ ... X^,Xr,„ ■ ■ ■ A'.JE) 
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(44) 

(45) 














by using (|^. The summation of the last four terms i.e. (41)+(42) is calculated as 


{l-2TL)i^{W\Xr, . Xr,_,\V), 


( 46 ) 


where (|^ is used. In total L cancels occur in (39) + (40) + (41) + (42). We have finished the 
proof. 
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